function plotrocks(filenumber,plottemphist)

    
N=readdata(0,1);
%     S=readdata(0,2);
%     T=readdata(0,3);
%     multigridlevels=readdata(0,4);
%     bound=readdata(0,5);
time=readdata(0,6);
%     dtnom=readdata(0,7);
dx=readdata(0,8);
dz=readdata(0,9);
mx=readdata(0,10);
mz=readdata(0,11);
%     mcs=readdata(0,12);
%     mct=readdata(0,13);
%     mtype=readdata(0,14);
%     mTemp=readdata(0,15);
%     mstressxx=readdata(0,16);
%     mstressxz=readdata(0,17);
%     mstrainxx=readdata(0,18);
%     mstrainxz=readdata(0,19);
%     mfinitestrain=readdata(0,20);
%     vx=readdata(0,21);
%     vz=readdata(0,22);
%     P=readdata(0,23);
Nsurfaces=readdata(0,24);
surfaces=readdata(0,25);

mcs=zeros(1,N); mct=zeros(1,N);
xdivision=1; xplotcellsize=sum(dx)/xdivision;
x=0; for s=0:xdivision-1 ind=find(mx>=x & mx<x+xplotcellsize); mcs(ind)=s; x=x+xplotcellsize; end
zdivision=150; zplotcellsize=sum(dz)/zdivision;
z=0; for t=0:zdivision-1 ind=find(mz>=z & mz<z+zplotcellsize); mct(ind)=t; z=z+zplotcellsize; end
Ia=int32(find((-1).^mcs.*(-1).^mct>0));
Ib=int32(find((-1).^mcs.*(-1).^mct<=0));
clear('mcs','mct');


for count=filenumber
    figure(1); clf;
    N=readdata(count,1);
    S=readdata(count,2);
    T=readdata(count,3);
    dx=readdata(count,8);
    dz=readdata(count,9);
    mx=readdata(count,10);
    mz=readdata(count,11);
    Nsurfaces=readdata(count,24);
    surfaces=readdata(count,25);

    x=zeros(1,S); z=zeros(1,T);
    for s=1:S-1 x(s+1)=x(s)+dx(s); end
    for t=1:T-1 z(t+1)=z(t)+dz(t); end

    minsurf=min(surfaces(:,Nsurfaces)); maxsurf=max(surfaces(:,Nsurfaces));

    FTmarkers=find(mz>=minsurf-0.5*min(dz) & mz<=maxsurf+0.5*min(dz));
    
    if plottemphist
        figure(2); clf;
        
        plot(x,surfaces(:,4),'r','linewidth',2); set(gca,'xlim',[0 2000e3]);
        hold on;
        for sur=4:1:Nsurfaces
            if mod(sur,5)==0 plot(x,surfaces(:,sur),'k','linewidth',0.5); end
        end
        plot(x,surfaces(:,Nsurfaces),'k','linewidth',2);
        
        plot(mx(FTmarkers),mz(FTmarkers),'.');
        set(gca,'ydir','reverse');
        %hold on;
        %plot(x,surfaces(:,Nsurfaces),'r'); set(gca,'ydir','reverse');
        picks=ginput;
        [Nhistpicks nothing]=size(picks); histpicks=zeros(1,Nhistpicks);
        for p=1:Nhistpicks
            [nothing histpicks(p)]=min((mx(FTmarkers)-picks(p,1)).^2+(100*(mz(FTmarkers)-picks(p,2))).^2);
        end
        plot(mx(FTmarkers(histpicks)),mz(FTmarkers(histpicks)),'r.'); drawnow;
        figure(1);
    end

    timehist=zeros(1,count); temphist=zeros(length(FTmarkers),count);
    for t=1:count
        t
        time=readdata(t,6);
        mTemp=readdata(t,15);
        timehist(t)=time/(60*60*24*365*1e6);
        temphist(:,t)=mTemp(FTmarkers);
    end
    timehist=-timehist+timehist(end);
    timehist=fliplr(timehist); temphist=fliplr(temphist);

    mFTage=zeros(1,length(FTmarkers));
    pdfAxis=linspace(0,20,1); lgeol=14.56; kinPar=4; pctPerTimeStep=1; annealModel=0; doProject=0; usedCf=0; kinParType=3; l0model=0; l0user=16;
    for e=1:length(FTmarkers)
        [cdf pdf oldestModelAge ftModelAge numPopulations]=ft(timehist,temphist(e,:),lgeol/l0user,kinPar,pctPerTimeStep,annealModel,doProject,usedCf,kinParType,l0model,l0user,pdfAxis);
        mFTage(e)=ftModelAge;
    end
    
    if plottemphist
        for p=1:Nhistpicks
            figure(10+p); clf;
            pdfAxis=linspace(0,20,100);
            [cdf pdf oldestModelAge ftModelAge numPopulations]=ft(timehist,temphist(histpicks(p),:)',lgeol/l0user,kinPar,pctPerTimeStep,annealModel,doProject,usedCf,kinParType,l0model,l0user,pdfAxis);
            subplot(2,1,1); plot(timehist,temphist(histpicks(p),:)'); title(['FT age: ' num2str(mFTage(histpicks(p))) ' Ma ']); xlabel('Time BP (Myr)'); ylabel('Temperature (^{\circ}C)');
            subplot(2,1,2); plot(pdfAxis,pdf); xlabel('FT Length (\mum)'); ylabel('Distribution');
        end
        figure(1);
    end
    
    xnew=linspace(min(x),max(x),5*S); surfnew=interp1(x,surfaces(:,Nsurfaces)',xnew);

    %FTprofile=griddata(mx(FTmarkers),mz(FTmarkers),mFTage,x,surfaces(:,Nsurfaces)');
    FTprofile=griddata(mx(FTmarkers),mz(FTmarkers),mFTage,xnew,surfnew);
    subplot(3,1,1);
    plot(xnew,FTprofile); set(gca,'xlim',[0 2000e3]); title([num2str(time/(60*60*24*365*1e6)) ' Ma ']);
    load ftprof_1; hold on; plot(x+907e3,ft_age_x,'*r');
    
    
    N=readdata(count,1);
    S=readdata(count,2);
    T=readdata(count,3);
%     multigridlevels=readdata(count,4);
%     bound=readdata(count,5);
    time=readdata(count,6);
%     dtnom=readdata(count,7);
    dx=readdata(count,8);
    dz=readdata(count,9);
    mx=readdata(count,10);
    mz=readdata(count,11);
%     mcs=readdata(count,12);
%     mct=readdata(count,13);
    mtype=readdata(count,14);
    %mTemp=readdata(count,15);
    %mstressxx=readdata(count,16);
    %mstressxz=readdata(count,17);
    %mstrainxx=readdata(count,18);
    %mstrainxz=readdata(count,19);
    %mfinitestrain=readdata(count,20);
%     vx=readdata(count,21);
%     vz=readdata(count,22);
     %P=readdata(count,23);
     Nsurfaces=readdata(count,24);
     surfaces=readdata(count,25);

    x=zeros(1,S); z=zeros(1,T); 
    for s=1:S-1 x(s+1)=x(s)+dx(s); end
    for t=1:T-1 z(t+1)=z(t)+dz(t); end


    
mTemp=readdata(count,15);
[grid,xgrid,zgrid]=bilingrid(mTemp,mx,mz,dx,dz,1,1); ind1=find(grid<1);
contvals=1200:25:1700; cont=contourc(xgrid,zgrid,grid',contvals);
clear('mTemp');


h1=subplot(3,1,3); mtype(find(mtype==6))=1; 
mtype(Ib)=mtype(Ib)+0.5;
[grid,xgrid,zgrid]=bilingrid(mtype,mx,mz,dx,dz,2,2); ind0=find(grid<0.75);
[Zgrid,Xgrid] = meshgrid(zgrid,xgrid); grid(ind0)=nan; pcolor(Xgrid,Zgrid,grid); shading flat; set(gca,'ydir','rev'); axis([0 2000000 0 710000]); caxis([0 5.5]); hold on;
hold on;
xlabel('Distance (m)'); ylabel('Depth (m)');


pos=1; allypoints=cell(1,length(contvals));
while pos<length(cont)
    step=cont(2,pos);
    i=find(contvals==cont(1,pos));
    xcont=cont(1,pos+1:pos+step); ycont=cont(2,pos+1:pos+step);
    plot(h1,xcont,ycont,'k','markersize',0.1);
    pos=pos+step+1;
    allypoints(i)={[cell2mat(allypoints(i)) ycont]};
end

subplot(3,1,2);
sealevel=0.5*(surfaces(1,4)+surfaces(S,4));
plot(x,-(surfaces(:,4)-sealevel),'r','linewidth',2); set(gca,'xlim',[0 2000e3]); set(gca,'ylim',[-5000 8000]);
hold on; 
for sur=4:1:Nsurfaces
    if mod(sur,5)==0 plot(x,-(surfaces(:,sur)-sealevel),'k','linewidth',0.5); end
end
plot(x,-(surfaces(:,Nsurfaces)-sealevel),'k','linewidth',2);
set(gca,'xticklabel','');
ylabel('Topography (m)');

print(gcf, '-djpeg100','-r600',['..\outfiles\rocksurf_ft' num2str(count)]);
 
end
end